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Abstract. The maximum-likelihood method for quantum estimation is reviewed 
and applied to the reconstruction of density matrix of spin and radiation as well 
as to the determination of several parameters of interest in quantum optics. 



1. Introduction 

Quantum estimation of states, observables and parameters is, from very basic 
principles, matter of statistical inference after sampling a population. Why is it 
so ? And what does this statement exactly mean ? The quantum description of a 
physical system is intrinsically a statistical one since given an unknown quantum 
state g : i) one cannot determine the quantum state from a single measurement 
(even joint measurement) performed on the system [[!]; ii) it is not possible to mea- 
sure different observables after copying the state preparation at disposal, since the 
linearity of quantum mechanics leads to the impossibility of cloning an arbitrary 
state [0. In a way, this means that any quantum estimation procedure necessary 
requires many identical preparations of the system under examination, on which 
one performs repeated measurements of an observable or of a set of observables. 

The full inference of the quantum state from feasible measurements is a hot 
topic of the last decade. In recent years, experiments have been performed to 
characterize the quantum state of different physical systems, such as single-mode 
radiation field [0] , a diatomic molecule [Q , a trapped ion Q , and an atomic beam 
Q. These results stimulated some efforts for efficient data-processing algorithms, 
in order to extract the maximum available information on the quantum state, since 
one always deals with finite ensembles [Q, and every detection scheme is affected 
by noise and imperfections. The need for an efficient data processing is even more 
stringent in cases where one is not interested in the complete characterization 
of the state, but only in some specific feature, as when one addresses the prob- 
lem of characterizing a device, rather than a quantum state, as, for example, the 
estimation of coupling constants, gain coefficients, or nonlinear susceptibilities. 

The most comprehensive quantum estimation procedure is quantum tomog- 
raphy ||. In quantum tomography the expectation value of an operator is ob- 
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tained by averaging a special function (usually termed sampling kernel or pattern 
function) over the experimental data of a sufficiently complete set of observables 
which is called a "quorum" . For example, in homodyne tomography of radiation 
the quorum observables are the quadratures of the e.m. field (for varying phase 
with respect to the local oscillator). The expectation value of a generic opera- 
tor is obtained by averaging the corresponding pattern function over data. The 
method is very general and efficient, however, in the averaging procedure, we have 
fluctuations which result in large statistical errors. 

In this paper we review the the maximum-likelihood (ML) principle approach 
to the quantum estimation problem. The ML idea is to find the quantum state, or 
the value of the parameters, that are most likely to generate the observed data. 
This idea can be quantified and implemented using the concept of the likelihood 
functional. Concerning the estimation of quantum state, in contrast to quantum 
tomography, the ML method estimates the state as a whole. As a result, a priori 
knowledge about properties of the density matrix can be incorporated from the 
very beginning, thus assuring positivity and normalization of matrix, with the re- 
sult of a substantial reduction of statistical errors ||. Regarding the estimation 
of specific parameters, we notice that in all the cases here analyzed the resulting 
estimators are efficient, unbiased and consistent, thus providing a statistically re- 
liable determination po|| . Moreover, by using the ML method only small samples 
of data are required for a precise determination. 

2. Maximum likelihood principle 

Here we briefly review the theory of the maximum-likelihood (ML) estimation of 
a single parameter. The generalization to several parameters, as for example the 
elements of the density matrix is, in principle, straightforward. The only point 
that should be carefully analyzed is the parameterization of the multidimensional 
quantity to be estimated. In the next section the specific case of the density matrix 
will be discussed. 

Let p(x\X) the probability density of a random variable x, conditioned to the 
value of the parameter A. The form of p is known, but the true value of the 
parameter A is unknown, and will be estimated from the result of a measurement 
of x. Let xi, X2, %n be a random sample of size N. The joint probability density 
of the independent random variable x\, X2, xn (the global probability of the 
sample) is given by 

C(xi,x 2 , ...,xn\X) = Ii^ =1 p(xk\X) , (1) 

and is called the likelihood function of the given data sample (hereafter we will 
suppress the dependence of C on the data). The maximum-likelihood estimator 
(MLE) of the parameter A is defined as the quantity A ML = A ML ({a;/ c }) that maxi- 
mizes C(X) for variations of A, namely A M l is given by the solution of the equations 
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Since the likelihood is positive the first equation is equivalent to OL/dX = where 



N 



L(A) = log£(A)=^logp(z fe |A) 



(3) 



fe=i 



is the so-called log-likelihood function. 

In order to obtain a measure for the confidence interval in the determination 
of A ML we consider the variance 



J[dx k p(x k \X) 



[A ML (K» - A] 



Upon defining the Fisher information 



F = I (Li- 
lt is easy to prove O] that 



d P (x\X) 



p(x\X) ' 



al^iNF)- 1 



(4) 



(5) 
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where N is the number of measurements. The inequality in Eq. (|6|) is known as the 
Cramer- Rao bound jfj) on the precision of ML estimation. Notice that this bound 
holds for any functional form of the probability distribution p(x\X), provided that 
the Fisher information exists VA and d\p(x\X) exists Wx. When an experiment 
has "good statistics" (i.e. a data sample large enough) the Cramer- Rao bound is 
saturated. 



3. Quantum state estimation 

In this section we review the the method for the maximum likelihood estimation 
of the quantum state, focusing attention to the cases of homodyne and spin tomo- 
graphies The physical situation we have in mind is an experiment consisting 
of N measurements performed on identically prepared copies of a given system. 
Quantum mechanically, each measurement is described by a positive operator- 
valued measure (POVM). The outcome of the ith measurement corresponds to 
the realization of a specific element of the POVM used in the corresponding run. 
We denote this element by IF. The likelihood is here a functional of the density 
matrix C(g) and is given by the product 

N 

£(g)=Y[^m), (7) 

i=l 

which represents the probability of the observed data. The unknown element of 
the above expression, which we want to infer from data, is the density matrix 
describing the measured ensemble. 
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We restrict ourselves to finite dimensional Hilbert spaces. In this case, it can 
be proved that C(g) is a concave function defined on a convex and closed set of 
density matrices. Therefore, its maximum is achieved either on a single isolated 
point, or on a convex subset of density matrices. In the first case we have a proper 
reconstruction scheme, namely the set of observables chosen for the measurement 
provides the complete characterization of the state under examinations. On the 
other hand, if the maximum is not unique, the set of chosen observables is insuf- 
ficient, namely it does not constitute a quorum. 

In order to optimize the procedure for the maximization of the likelihood func- 
tion we introduce a specific parameterization of the set of density matrices. A 
given density matrix can be written in the form 

g = ftf, (8) 

which automatically guarantees that g is positive and Hermitian for any com- 
plex lower triangular matrix T, with real elements on the diagonal. For an M- 
dimensional Hilbert space, the number of real parameters in the matrix T is 
M + 2 M (M— 1) /2 = M 2 , which equals the number of independent real parameters 
for a Hermitian matrix. This confirms that our parameterization is minimal, up 
to the unit trace condition. 

In numerical calculations, it is convenient to replace the likelihood functional 
by its natural logarithm, which does not change the location of the maximum. 
Thus the function subjected to numerical maximization is given by 

N 

L(f) = 53 In Tr(f + fni) - ATr(f t T) , (9) 

i=l 

where A is a Lagrange multiplier accounting for normalization of g that equals the 
total number of measurements N []l3[ . Using this formulation, the maximization 
problem can be solved by standard numerical procedures for searching the max- 
imum over the M 2 real parameters of the matrix T ||. The examples presented 
below use the downhill simplex method Jl4| . 

Our first example is the application of the ML estimation in quantum homo- 
dyne tomography of a single- mode radiation field p5[ , which is so far the most 
successful method in measuring nonclassical states of light |5|,|l6| . The experimental 
apparatus used in this technique is the homodyne detector. The realistic, imperfect 
homodyne measurement is described by the positive operator-valued measure 

H(x; <p) = - 1 exp ( J X ~ Vvi V ) 2 \ {1Q) 

vM 1 -v) V 1 -v J 

where r\ is the detector efficiency, and x v = (a* e ttp + ae~ lip )/2 is the quadrature 
operator ([a, a'] = 1), depending on the externally adjustable local oscillator (LO) 
phase ip. 

After repeating the measurement N times, we obtain a set of pairs (x^ipi) 
consisting of the outcome Xi and the LO phase ipi for the iih run, where i = 
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1, . . . , N. The log-likelihood functional is given by Eq. (|) with ft; = H(xi] <{*). Of 
course, for a light mode it is necessary to truncate the Hilbcrt space to a finite 
dimensional basis. We shall assume that the highest Fock state has M — 1 photons, 
i.e. that the dimension of the truncated Hilbert space is M. For the expectation 
Tr[ftf H(x; <p)\ it is necessary to use an expression which is explicitly positive, 
in order to protect the algorithm against occurrence of small negative numerical 
arguments of the logarithm function. A simple derivation yields 



M-l k 
k=0 j=0 



k-j 

J2(k\T\n + j)B n+J , n (n\x)e m * 

n=0 



(11) 



where B n+J . n = [("+ J >"(1 - and (n\x) = H n {x) exp(-a?/2)/y/V>nW 

are eigenstates of the harmonic oscillator in the position representation — H n (x) 
being the nth Hermite polynomial. 

We have applied the ML technique to reconstruct the density matrix in the Fock 
basis from Monte Carlo simulated homodyne statistics § . Fig. [l] depicts the matrix 
elements of the density operator as obtained for a coherent state and a squeezed 
vacuum. Remarkably, only 50000 homodyne data have been used for quantum 
efficiency at photodctcctors r\ = 80%. We notice that the ML method is affected 
by much smaller statistical errors than conventional tomography. As a comparison 
one could see that the same precision of the reconstructions in Fig. [l] could be 
achieved using 10 7 -10 8 data samples in conventional tomography of Ref. |p!5[ . On 
the other hand, in order to find numerically the ML estimate we need to set a 
priori the cut-off parameter for the photon number, and its value is limited by 
increasing computation time. 



Po, 





Figure 1. Reconstruction of the density matrix of a single-mode radiation field by the ML 
method. The plot shows the matrix elements of a coherent state (left) with {fit a) = 1 photon, and 
for a squeezed vacuum (right) with (a' a) = 0.5 photon. A sample of 50000 simulated homodyne 
data for quantum efficiency r] = 80% has been used. 

We mention that ML estimation can also be applied to the reconstruction of 
the quantum state a two- mode field Q, along with the multi-mode tomographic 
technique with a single LO Jl7[ | . 

Finally, we apply the ML procedure for reconstructing the density matrix of 
spin systems. For example, let us consider TV repeated preparations of a pair of 
spin-1/2 particles. The particles are shared by two parties. In each run, the parties 
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select randomly and independently from each other a direction along which they 
perform spin measurement. The obtained result is described by the joint projection 
operator over spin coherent states 

£ = |fif,Qf)<Qf,Qf|, (12) 

where £lf and fif are the vectors on the Bloch sphere corresponding to the out- 
comes of the iih run, and the indices A and B refer to the two particles. As 
in the previous examples, it is convenient to use an expression for the quantum 
expectation value Tr(T'T^-i) which is explicitly positive. The suitable form is 

Tr(ftf^) = ^|< M |f|nf ) of)| 2 , 

where is an orthonormal basis in the Hilbcrt space of the two particles. The 
result of a simulated experiment with only 500 data for the reconstruction of the 
density matrix of the singlet state is shown in Fig. 0. 




Figure 2. Reconstruction of the density matrix of a pair of spin- 1/2 particles in the singlet state 
by ML method. The matrix elements has been obtained by a sample of 500 simulated data. 



4. Parameters estimation in quantum optics 

Here we focus our attention on the determination of specific parameters which are 
relevant in quantum optics, and analyze their ML estimation procedure in some 
details. In the next two subsections we consider the estimation of the parameters 
of a Gaussian state and the estimation of the quantum efficiency of both linear 
and avalanche photodetectors. The reader may found more details in Ref. |To|| . 

4.1. GAUSSIAN STATE ESTIMATION 

In this section we apply the ML method to estimate the quantum state of a 
single-mode radiation field that is characterized by a Gaussian Wigner function. 
Such kind of states comprises the wide class of coherent, squeezed and thermal 
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states, namely most of the states effectively produced in an optical laboratory. We 
consider the Wigner function of the form 



2A 



2 



W{x, y) = exp <^ -2A 



,2 



~(x- a) 2 + n(y - bf 

K 



(13) 



and we apply the ML technique starting from homodyne detection to estimate 
the four real parameters A, k, a and b. The four parameters are connected to the 
number of thermal, squeezing and coherent-signal photons in the quantum state 
as follows 

1 (JL _ A - 1 + k2 _ I 

2 I A 2 Usq ~ 4 K 2 



n *h = n ( 72 _ 1 ) n sq = — ^ n coh = a 2 + b 2 . (14) 



In fact, the quantum state corresponding to the Wigner function in Eq. (|T^) writes 

§ = D(riS(r)-±- (-^-rY " S\r)DH P ) , (15) 
n th + I \n t h + i- J 

with r — \ log k and /i = a + ib, and where S(r) = exp[r(a 2 — a^ 2 )/2] and D(/i) = 
exp(fj,cv — n*a) denote the squeezing and displacement operators, respectively. 

We consider repeated preparations of a Gaussian state, on which we perform 
homodyne measurements at different phases <j) with respect to the local oscillator. 
The homodyne distribution is given, for unit quantum efficiency of photodetectors, 
by the Gaussian M 



exp\-2A 2 K^- / , 2 Vn . (16) 

sin ffl) n 1 cos^ <b + sin <b 



For non-unit quantum efficiency the ideal distribution of Eq. (|lq ) is replaced by a 
convolution with a Gaussian of variance (l-J7)/(4»7). From Eqs. (|) and ( |l6|) one 
easily evaluates the log-likelihood function for a set of N homodyne outcomes Xi 
at random phases 4>i a s follows 

L = y I l og 2A ' K 2 - 2A2K [^-"CQs0-6sin0)] 2 

2 7t(k 2 cos 2 0i + sin (pi) k 2 cos 2 </>i + sin fa 

The ML estimators A M l , Kml , Oml and 6 ML are found upon maximizing Eq. ( |l7| ) 
versus A, k, a and b. 

In order to check the reliability of the state reconstruction we performed a 
set of Monte Carlo simulated experiments, starting from homodyne measurements 
with quantum efficiencies in the range n — 70 — 90%. For states with average 
photons in the range n t h < 3, n co h < 5, and n sq < 3 and for data samples of size 
of the order N — 10 4 — 10 5 we always found a reconstructed state very close to 
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the theoretical one. [] As a matter of fact, the quality of the state reconstruction is 
good enough that other physical quantities that are theoretically evaluated from 
the experimental values of A ML ,K ML ,a ML and b UL are inferred very precisely. For 
example, we evaluated the photon number probability of a squeezed thermal state, 
which is given by the integral 

/ i -i \ y 3 *^ to, "«/>,*)-i] w MR x 

W<A») = J c^ nth , K y^ ■ (18) 

with C(4>, nth, w) = (nth + sin + kcos 2 <j>) + A. The comparison of the 

theoretical and the experimental results for a state with nth = 0.1 and n sq = 3 is 
reported in Fig. ||. The statistical error of the reconstructed number probability 
affects the third decimal digit, and is not visible on the scale of the plot. 
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Figure 3. The photon distribution of a squeezed-thermal state as reconstructed by the ML 
estimation of the corresponding Wigner function (left histogram for the theoretical values, right 
histogram for the reconstructed values). Number of data samples N = 50000, quantum efficiency 
r) = 80%, number of thermal photons nth = 0.1, number of squeezing photons n sq = 3. The 
statistical error affects the third decimal digit, and it is not visible in the scale of the plot. 



As a further development, we mention that the ML estimation of Gaussian 
Wigner functions also provides a technique to estimate the coupling constants of 
quadratic Hamiltonians of the form 

H = aa + a*a) + ipa^a + ^£a 2 + -£*a t2 . (19) 



Hamiltonians like that in Eq. (19) describes the interaction of light modes in active 
optical medium characterized by a second order susceptibility tensor. Actually, the 



1 As a global measure of the goodness of the reconstruction one should consider the normalized 
overlap O between the theoretical and the estimated state 

_ Tr \eSm\ 



VTrMTrHj 



which is unit O = 1 if and only if q = q M l ■ 
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unitary evolution operator U — e~ lHt preserves the Gaussian form of an input 
state with Gaussian Wigner function, and therefore one can use a Gaussian state 
to probe and characterize an optical device. 

4.2. ABSOLUTE ESTIMATION OF THE QUANTUM EFFICIENCY 

The operation of a photodetector is, in principle, very simple: each photon ionizes 
an atom, and the resulting charge is amplified to produce a measurable pulse. In 
practice, however, available photodetectors are usually characterized by a quan- 
tum efficiency lower than unity, which means that only a fraction of the incoming 
photons lead to an electric pulse, and ultimately to a " count" . We may distinguish 
two main classes of photodetectors. In the first we have detectors where the result- 
ing current proportional to the incoming photon flux: in this case we have a linear 
photodetector. For example, this is the case of the high flux photodetectors used in 
homodyne detection. In the second class we have photodetectors operating at very 
low intensities, which resort to avalanche process in order to transform a single 
ionization event into a recordable pulse. This implies that one cannot discriminate 
between a single photon or many photons as the outcomes from such detectors 
are either a " click" , corresponding to any number of photons, or " nothing" which 
means that no photons have been revealed. 

Conventional characterization of photodetectors resorts to prepare a reference 
state with known intensity, and then measuring which fraction of the signal is ac- 
tually revealed. This unavoidably leads to rather poor performances when applied 
in the relevant regime of quantum signals. Detection losses, in facts, distorce the 
whole probability distribution of the quantity being measured, not only the aver- 
age value. Moreover, we need the accurate knowledge of the quantum state of the 
reference signal. In the following, we apply the ML principle to the absolute esti- 
mation of the quantum efficiency of both linear and avalanche photodetectors. We 
show that, along with the reliable characterization of quantum signals, ML method 
is an effective and statistically efficient tool for characterizing the response of a 
photodetector to low-intensity and/or nonclassical states. 

Let us first study the case of linear photodetectors. As a reference state we 
consider a squeezed-coherent state, measured by homodyne detection. The effect 
of non-unit quantum efficiency r\ on the probability distribution of homodyne 
detection is twofold. We have both a rescaling of the mean value and a broadening 
of the distribution. For a squeezed state \xo,r) — D(xo)S(r)\0) with the direction 
of squeezing parallel to the signal phase and to the phase of the homodyne detection 
(without loss of generality we set this phase equal to zero and xq , r > ) we have 




1 



(x - 1]X ) 

2A 2 



2 1 



1 



p v {x) 



exp 




(20) 



The total number of photons of the state is given by n = Xq + sinh 2 r, whereas the 
squeezing fraction is defined as 7 = sinh 2 r/n. Apart from an irrelevant constant, 
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the log-likclihood function can be written as 

- L(r?) = log A 2 + i (x 2 + rjxl - 2r]x x 
The resulting MLE is thus given by 



(21) 



= 1 




Mx 2 



+ (1 + e- 2r ){x - 2x + x e- 2r )x 



A set of Monte Carlo simulated experiments confirmed that the Cramer-Rao 
bound is attained. The performances of the ML estimation can be compared to 
the "naive" estimation based only on the measurement of the mean value, i.e. 
f]A\ — x/xq. We expect the naive method to be less efficient, since the quantum 
efficiency not only rescales the mean value, but also spreads the variance of the ho- 
modyne distribution in Eq. (|2C|). In Fig. [|, on the basis of a Monte Carlo simulated 
experiment, we compare the ML and the average-value methods in estimating the 
quantum efficiency through homodyne detection on a squeezed state. The advan- 
tages of ML method are apparent, especially for the estimation of low values of r\. 
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Figure 4- Quantum efficiency of linear photodetectors: ML estimation through homodyne 
detection on a squeezed state. The plots report the ratio between the estimated value of the 
quantum efficiency and the true value, as a function of the true value. On the left the results 
from maximum-likelihood method; on the right by the "naive" average- value method. Homodyne 
sample: 2500 data. Reference state: a squeezed state with mean number of photons n = 1 and 
squeezing fraction 7 = 99% (nearly a squeezed vacuum). 



Let us now consider avalanche photodetectors, which perform the ON/OFF 
measurement described by the two-value POM 

00 

II OFF = - V) P \P) (P| Hon = I - IIoff , (22) 

where I denotes the identity operator. With avalanche photodetectors we have 
only two possible outcomes: "click" or "no clicks" which we denote by "1" and "0" 
respectively. The log-likclihood function is given by 



L(n) = {N- N c ) logP (?7) + N c log[l - P fa)] , 



(23) 
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where Po(r)) — Tr[^n FF] is the probability of having no clicks for the reference 
state described by the density matrix g, N is the total number of measurements, 
and N c is the number of events leading to a click. The maximum of L(rj), i.e. the 
MLE for the quantum efficiency, satisfies the equation 

Po(% L ) = l-§, (24) 

whose solution, of course, depends on the choice of the reference state. The optimal 
choice would be using single-photon states as a reference. In this case, we have the 
trivial result i] ML — N c /N. However, single-photon state are not easy to prepare 
and generally one would like to test r\ for coherent pulses |a). In this case, we have 
Po(v) — ex P(—\ a \ 2 v) an d 

%L = -^log(l-§) . (25) 

The Fisher information is given by 

F= (Sr) ^ + (Sr) ^ = p (i-po) (w) ' (26) 

and therefore, for a weak coherent reference one has 

e vH 2 _i | a |2 i V ' 



5. Summary and conclusions 

In this paper we reviewed the application of the maximum likelihood principle 
to the reconstruction of the density matrix of a generic quantum system ||, as 



well as to the estimation of relevant parameters in quantum optics 10 . In all 
cases, the resulting reconstruction algorithm is statistically efficient, and provides 
the reliable estimation of the quantity of interest using much smaller data samples 
compared to conventional methods. In particular, the ML estimation of the density 
matrix allows one to incorporate the natural physical constraints we have on the 
quantum state, thus leading to a substantial reduction of statistical fluctuations. 
For quantum-optical parameters, the ML approach provides efficient estimation 
schemes based on feasible measurements like homodyne detection. The resulting 
procedures lead to substantial improvement over conventional methods and are of 
technological interest. 
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